// Revision 3: Simplify framework.
#include "OpenCVSample0.h"

double x00[30];
double alpha=1;
double w[7];
double rate[7],brate[7];

int itr=0;

#define FC 31

#define SVM_CSV "..\\resources\\svm.csv"
#define TRAIN_DATA "..\\resources\\trdata.xml"
#define MASK_JPG "..\\resources\\mask.jpg"
#define MASK3_JPG "..\\resources\\mask3.jpg"
#define IMAGE_H "..\\resources\\r34_h.jpg"
#define IMAGE_M "..\\resources\\r34_m.jpg"
#define IMAGE_D "..\\resources\\r34_d.jpg"

static void read_csv(const string& filename, vector<Mat>& images, vector<float>& labels, char separator = ';') {
    std::ifstream file(filename.c_str(), ifstream::in);
    if (!file) {
        string error_message = "No valid input file was given, please check the given filename.";
        CV_Error(CV_StsBadArg, error_message);
    }
    string line, path, classlabel;
    while (getline(file, line)) {
        stringstream liness(line);
        getline(liness, path, separator);
        getline(liness, classlabel);
        if(!path.empty() && !classlabel.empty()) {
            images.push_back(imread(path,CV_LOAD_IMAGE_COLOR));
            labels.push_back(atof(classlabel.c_str()));
	    //printf("PUSH:%d\n",atoi(classlabel.c_str()));
        }
    }
}
void  function1_fvec(const real_1d_array &x, real_1d_array &fi, void *ptr)
{
    CvSVM SVM;
    SVM.load(TRAIN_DATA);
    Mat ProjDest(1,FC,CV_32FC1);
    double move[4];
    double w[4];
    double tmove=0;
    //printf("====================ITR %d==============\n",itr++);
    double d=0.0001;

    
    //weighted move
    double r[7];
   
    for(int j=0;j<7;j++)
	r[j]=1;
    r[1]=3;
    r[4]=10;r[5]=30;r[6]=10;
    for(int j=0;j<7;j++)
    {
	tmove+=r[j]*pow(x[j]-rate[j],2);
    }
   
    for(int j=0; j<7; j++)
        ProjDest.at<float>(j)=x00[j]*(1+x[0]);
    for(int j=7; j<14; j++)
        ProjDest.at<float>(j)=x00[j]*(1+x[1]);
    for(int j=14; j<21; j++)
        ProjDest.at<float>(j)=x00[j]*(1+x[2]);
    for(int j=21; j<28; j++)
        ProjDest.at<float>(j)=x00[j]*(1+x[3]);
    ProjDest.at<float>(28)=x00[28]*(1+x[4]);
    ProjDest.at<float>(29)=x00[29]*(1+x[5]);
    ProjDest.at<float>(30)=x00[30]*(1+x[6]);
    
    //standard derivation
    double mean=0,stdev=0;
    for(int j=0;j<4;j++)
    {
	mean+=x[j];
	stdev+=pow(x[j],2);
    }
    mean=mean/4;
    stdev=stdev-mean*mean;
    for(int j=0;j<7;j++)
    {
	//printf("x[%d]->%lf ",j,x[j]);
    }
    //printf("\n");
    
    fi[0]=-SVM.predict(ProjDest)+alpha*sqrt(tmove)+5;
    //fi[0]=10/(-0.7*SVM.predict(ProjDest)+move*0.3);
    //printf("%lf,SVM=%lf,move=%lf,stdev=%lf\n",fi[0],SVM.predict(ProjDest),tmove,stdev);
}

double calcEnt(Mat &M,Mat mask)
{
    vector<Mat> mv;
    split(M,mv);
    long hist[256];
    for(int j=0;j<256;j++)
	hist[j]=0;
    long pixCount=0;
    for(int x=0;x<mv[0].rows;x++)
    {
	for(int y=0;y<mv[0].cols;y++)
	{
	    if(mask.at<Vec3b>(x,y)[0])
	    {
		int r=mv[0].at<uchar>(x,y);
		//printf("%d %d\n",r,hist[r]);
		pixCount+=1;
		hist[r]+=1;
	    }
	}
    }
    double ent=0;
    printf("%ld\n",pixCount);
    for(int j=0;j<256;j++)
    {
	double r=hist[j]/(double)pixCount;
	//printf("%lf\n",r);
	if(r!=0)
	    ent-=r*log(r)/log(2);
    }
    return ent;
}

Mat srcb,src_hb,src_mb;
Mat mask,maskR;

void scrollbar_onchange(int count,void *data)
{
	Mat src,src_h,src_m,dest;
	srcb.copyTo(src);
	src_hb.copyTo(src_h);
	src_mb.copyTo(src_m);
	//cout << count <<endl;
	double finalrate=(double)(count)/30.0-1;
	for(int i=0;i<7;i++)
		rate[i]=brate[i]*(1+finalrate);
	clock_t totaltime=clock();
	#pragma omp parallel
	{

		/*while(currEnt>destEnt&&itercount<35)
		{
		GaussianBlur( src_h, dest, Size( 2*itercount+1, 2*itercount+1 ), 0, 0 );
		currEnt=calcEnt(dest,mask);
		itercount+=1;
	   
		}
		printf("currEnt=%lf,K=%d\n",currEnt,2*itercount-1);
		if(itercount!=0)
		GaussianBlur(src_h, src_h, Size( 2*itercount-1, 2*itercount-1 ), 0, 0 );*/
		#pragma omp sections 
		{
			#pragma omp section
			{
				clock_t begin_time = clock();
	
				vector <Mat> mv;
				int size=20*fabs(rate[0]);
				//gpu::GpuMat gpuSrcH(src_h);
				//gpu::GaussianBlur(gpuSrcH,gpuSrcH,Size( 31, 31 ), 0, 0 );
				GaussianBlur(src_h, src_h, Size( 2*size+1, 2*size+1 ), 0, 0 );
				//std::cout << "Elapsed time / GaussianBlur: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
		
				//gpuSrcH.download(src_h);
				begin_time=clock();
				cvtColor(src_h,src_h,CV_BGR2HLS);
				split(src_h,mv);
				mv[2].convertTo(mv[2],-1,rate[1]+1,0);
				merge(mv,src_h);
				cvtColor(src_h,src_h,CV_HLS2BGR);
				//std::cout << "Elapsed time / TestCopy : "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
			}
			#pragma omp section
			{
				clock_t begin_time=clock();
				//printf("===full aFter===\n");
				// full image remove
				{
					Mat mask2(src_m.rows,src_m.cols,CV_8UC1);
					cvtColor(mask,mask2,CV_BGR2GRAY);
					Scalar stdev,mean;
					Mat dest(src_m.rows,src_m.cols,CV_8UC1);
					cvtColor(src_m,dest,CV_BGR2GRAY);
					meanStdDev(dest,mean,stdev,mask2);    
					double pixThres=mean.val[0]-2*stdev.val[0];
					double desThres=mean.val[0]-stdev.val[0];
					double pixMean[3];
					double pixStdev[3];
					meanStdDev(dest,mean,stdev,mask2); 
					//printf("Mean:%lf Stdev:%lf,Thres=%lf\n",mean.val[0],stdev.val[0],pixThres);
	
					float range[] = { 0, (int) pixThres } ;
					int histSize = (int) pixThres;
					const float* histRange = { range };
					bool uniform = true; bool accumulate = false;

					Mat hist;
					meanStdDev(mask2,mean,stdev); 
					double wcount=mean.val[0]/255.0*(src_m.rows*src_m.cols);
	
					calcHist( &dest, 1, 0, mask2, hist, 1, &histSize, &histRange, uniform, accumulate );
					//printf("%d\n",hist.rows);
					double aveG=0;
					int xt=0;
					double tG=0;
					for(int q=0;q<hist.rows;q++)
						tG+=hist.at<float>(q)/wcount;
	    
					for(int q=0;q<hist.rows;q++)
					{
						aveG+=hist.at<float>(q)/wcount;
						//printf("q:%d,%lf\n",q,aveG);
						if(aveG>=(1-fabs(rate[3]))*tG)
						{
							//printf("DESTRATE:%lf,q=%d\n",1-fabs(rate[3]),q);
							xt=q;
							break;
						}
		    
					}
					//std::cout << "Elapsed time / Histogram: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
    
					begin_time=clock();
					//Mat dest2;
					//dest.copyTo(dest2);
					#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						uchar *t;
						Vec3b *s;
						t=dest.ptr<uchar>(x);
						s=mask2.ptr<Vec3b>(x);
						Vec3b *w;
						w=src_m.ptr<Vec3b>(x);
						for(int y=0;y<src_m.cols;y++)
						{
							//printf("%d %d\n",x,y);
							if(t[y]<xt && s[y][0]!=0)
							{
								//t[y]=t[y]*pixThres/xt;
								for(int v=0;v<3;v++)
									w[y][v]=w[y][v]*pixThres/xt;
								//printf("Yeah: %lf\n",(double)pixThres/(double)xt);
							}
							if(t[y]>=xt && t[y]<pixThres && s[y][0]!=0)
							{
								//t[y]=pixThres;
								double dg=0.2126*w[y][0]+0.7152*w[y][1]+0.0722*w[y][2];
								for(int v=0;v<3;v++)
									w[y][v]=w[y][v]*pixThres/dg;
							}
						}
					}
					//getchar();
					/*#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						for(int y=0;y<src_m.cols;y++)
						{
							for(int t=0;t<3;t++)
							{
								if(dest2.at<uchar>(x,y)!=0)
									src_m.at<Vec3b>(x,y)[t]=src_m.at<Vec3b>(x,y)[t]*dest.at<uchar>(x,y)/dest2.at<uchar>(x,y);
								else
									src_m.at<Vec3b>(x,y)[t]=0;
							}
						}
					}*/
					//std::cout << "Elapsed time / frecks: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
    
					begin_time=clock();
					//dest.copyTo(dest2);
					meanStdDev(dest,mean,stdev,mask2); 
					double meanval=mean.val[0];
					#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						Vec3b *w;
						w=src_m.ptr<Vec3b>(x);
						uchar *v;
						v=dest.ptr<uchar>(x);
						for(int y=0;y<src_m.cols;y++)
						{
							double t=meanval-(meanval-v[y])*(1+rate[2]);
							//double dg=0.2126*w[y][0]+0.7152*w[y][1]+0.0722*w[y][2];
							for(int j=0;j<3;j++)
								w[y][j]=w[y][j]*t/v[y];

						}
					}
					/*#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						for(int y=0;y<src_m.cols;y++)
						{
							for(int t=0;t<3;t++)
							{
								if(dest2.at<uchar>(x,y)!=0)
									src_m.at<Vec3b>(x,y)[t]=src_m.at<Vec3b>(x,y)[t]*dest.at<uchar>(x,y)/dest2.at<uchar>(x,y);
								else
									src_m.at<Vec3b>(x,y)[t]=0;
							}
						}
					}*/
					//std::cout << "Elapsed time / stddev: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
				}
			}
		}
	 }
		//imshow("src_m",src_m);
		//waitKey(0);
		clock_t begin_time=clock();
		//printf("Maybe Here?\n");
		Mat dif=src_h;
		Mat ddif=src_m;
		Mat dst(dif.rows,dif.cols,CV_8UC3);
		double smax=393.088;
		double smin=-0.969;
		double dmax=97.104;
		double dmin=-76.391;
		double ddmax=90;
		double ddmin=-90;
		/*Mat doubleh(src_h.rows,src_h.cols,CV_32FC3);
		Mat doublem(src_h.rows,src_h.cols,CV_32FC3);
		Mat doubled(src_h.rows,src_h.cols,CV_32FC3);
		clock_t testtime=clock();
		#pragma omp parallel 
		{
			#pragma omp sections
			{
				#pragma omp section
				{
					clock_t sectime=clock();
					src_h.convertTo(doubleh,CV_32FC3);
					doubleh.convertTo(doubleh,-1,(dmax-dmin)/255,dmin);
					std::cout << "Elapsed time / total: "  << float( clock () - sectime ) /  CLOCKS_PER_SEC << endl;	
				}
				#pragma omp section
				{
					clock_t sectime=clock();
					src_m.convertTo(doublem,CV_32FC3);
					doublem.convertTo(doublem,-1,(ddmax-ddmin)/255,ddmin);
					std::cout << "Elapsed time / total: "  << float( clock () - sectime ) /  CLOCKS_PER_SEC << endl;	
				}
				#pragma omp section
				{
					clock_t sectime=clock();
					src.convertTo(doubled,CV_32FC3);
					doubled.convertTo(doubled,-1,(smax-smin)/255,smin);
					std::cout << "Elapsed time / total: "  << float( clock () - sectime ) /  CLOCKS_PER_SEC << endl;	
				}
			}
		}
		std::cout << "Elapsed time / total: "  << float( clock () - testtime ) /  CLOCKS_PER_SEC << endl;
		clock_t sectime=clock();
		add(doubleh,doublem,doubleh,Mat(),-1);
		std::cout << "Elapsed time / ADD: "  << float( clock () - testtime ) /  CLOCKS_PER_SEC << endl;
		add(doubleh,doubled,doubled,Mat(),-1);
			doubled.copyTo(dst);
		*/
		Vec3b *p,*q,*a,*b;
		#pragma omp parallel for num_threads(8) private(p,q,a,b)
		for(int i=0;i<src.rows;i++)
		{
		  for(int j=0;j<src.cols;j++)
		  {
			  p = src.ptr<Vec3b>(i);
			  q = dif.ptr<Vec3b>(i);
			  a = ddif.ptr<Vec3b>(i);
			  b = dst.ptr<Vec3b>(i);
			  for(int t=0;t<3;t++)
			  {
				double l=p[j][t];
				double r=q[j][t];
				double m=a[j][t];
				l=l/255*(smax-smin)+smin;
				r=r/255*(dmax-dmin)+dmin;
				m=m/255*(ddmax-ddmin)+ddmin;
				if(l+r+m>255)
					b[j][t]=255;
				else
				{
					if(l+r+m<0)
						b[j][t]=0;
					else
						b[j][t]=l+r+m;
				}
			  }
		   }
		}
		//std::cout << "Elapsed time / trimerge: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;
    
    //system("./trimerge ./r34_d.jpg ./final_h.jpg ./final_m.jpg 393.088 -0.969 97.104 -76.391 90 -90");
	begin_time=clock();
    Mat corrSrc=dst;
    cvtColor(corrSrc,corrSrc,CV_BGR2HLS);
	/*vector <Mat> mv;
	split(corrSrc,mv);
	#pragma omp parallel for
	for(int i=0;i<3;i++)
		mv[i].convertTo(mv[i],-1,1+rate[4+i],0);
	merge(mv,corrSrc);*/
	
	#pragma omp parallel for num_threads(8) private(p,q,a,b)
    for(int x=0;x<corrSrc.rows;x++)
    {
		p = corrSrc.ptr<Vec3b>(x);
		q = src.ptr<Vec3b>(x);
		for(int y=0;y<corrSrc.cols;y++)
		{
			if(maskR.at<Vec3b>(x,y)[0])
			{
				for(int j=0;j<3;j++)
				{
					int r=p[y][j];
					if(r*(1+rate[j+4])<255)
						p[y][j]=r*(1+rate[j+4]);
					else
						p[y][j]=255;
				}
			}
		}
    }
    cvtColor(corrSrc,corrSrc,CV_HLS2BGR);
	//std::cout << "Elapsed time / other: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
	//std::cout << "Elapsed time / total: "  << float( clock () - totaltime ) /  CLOCKS_PER_SEC << endl;
	imshow("Result",corrSrc);
}



Mat main_process( double alphaIn, double argv6, char* imagePath_h, char* imagePath_m, char* imagePath_d)
{
	//std::cout << "It works!" << endl;
    alpha=alphaIn;
    ifstream in(SVM_CSV);
    string line, field;

    vector< vector<double> > array;  // the 2D array
    vector<double> v;                // array of values for one line only

    while ( getline(in,line) )    // get next line in file
    {
        v.clear();
        stringstream ss(line);

        while (getline(ss,field,','))  // break line into comma delimitted fields
        {
            v.push_back(atof(field.c_str()));  // add each field to the 1D array
        }

        array.push_back(v);  // add the 1D array to the 2D array
    }
   
    // TODO:
    // 1 -> simplify framework
    // 2 -> gradient calculation
    // 3 -> optimization rebuild
    // 4 -> mask support put to r4
    // 5 -> change HSV value.
    // Color correction will not be considered in this version
    // This version will also be a playground for all new ideas.
    // El psy congroo.
    
    
    double features [34][FC];
    double labels [34];
    
    Mat trainingDataMat(34, FC, CV_32FC1); 
    Mat labelsMat(34, 1, CV_32FC1); 
    for(int i=0;i<34;i++)
    {
      for(int j=0;j<FC;j++)
      {
			trainingDataMat.at<float>(i,j)=array[i+1][j];
			//printf("%lf ",trainingDataMat.at<float>(i,j));
      }
      labelsMat.at<float>(i)=array[i+1][FC];
      //printf("%lf \n",labelsMat.at<float>(i));
    }
        
    // training 
    CvSVMParams params;
    
    params.svm_type    = CvSVM::EPS_SVR;
    params.kernel_type = CvSVM::RBF;
    params.term_crit   = cvTermCriteria(CV_TERMCRIT_ITER, 10000, 1e-7);
    params.nu     = 0;
    params.C      = 62.50;
    params.gamma  = 0.02250;
    params.p      = 0.01;

    // Train the SVM
    CvSVM SVM;
    //printf( "\nTraining the SVM (in progress) ..... ");
                                
                // train using auto training parameter grid search if it is available 
                // (i.e. OpenCV 2.0) with 10 fold cross valdiation
                // N.B. this does not search kernel choice 
   
    //SVM.train_auto(trainingDataMat, labelsMat, Mat(), Mat(), params);
    SVM.train(trainingDataMat, labelsMat, Mat(), Mat(), params);
    params = SVM.get_params();
    //printf( "\nUsing optimal parameters degree %f, gamma %f, ceof0 %f\n\t C %f, nu %f, p %f\n Training ..\n", params.degree, params.gamma, params.coef0, params.C, params.nu, params.p);
    SVM.save(TRAIN_DATA);
    //SVM.load("trdata_good.xml");
    
    Mat sampleMat(1,FC,CV_32FC1);
    Mat predictMat(1,FC,CV_32FC1);
    for(int i=0;i<34;i++)
    {
      for(int j=0;j<FC;j++)
      {
		sampleMat.at<float>(j)=array[i+1][j];
      }
      double response=SVM.predict(sampleMat);
      //printf("response=%lf,actual is %lf\n",response,array[i+1][128]);
     // printf("%lf %lf\n",array[i+1][FC],response);
    }
    
    // priori knowledge acquisation
    for(int j=0;j<FC;j++)
    {
	sampleMat.at<float>(j)=array[34][j];
    }
    double dxmax=0;
    double xmax=0;
    for(double x=0;x>-1;x=x-0.001)
    {
	for(int j=0;j<FC;j++)
	{
	    //printf("%lf ",sampleMat.at<float>(j));
	}
	//printf("\n");
	for(int t=0;t<7;t++)
	{
	    sampleMat.at<float>(t+14)=array[34][t+14]*(1.001+x);
	}
	double response=SVM.predict(sampleMat);
	for(int t=0;t<7;t++)
	{
	    sampleMat.at<float>(t+14)=array[34][t+14]*(1+x);
	}
	double response2=SVM.predict(sampleMat);
	double dx=(response2-response)/0.001;
	if(dx>dxmax)
	{
	    dxmax=dx;
	    xmax=x;
	}
	//printf("%lf\t%lf\t%lf\n",x,dx,response);
	
    }
    
    //printf("xmax=%lf,dx=%lf\n",xmax,dxmax);
    for(int j=0;j<FC;j++)
    {
	sampleMat.at<float>(j)=array[34][j];
    }
    double dymax=-1,ymax=-1;
    for(double x=0;x>-1;x=x-0.001)
    {
	for(int j=0;j<FC;j++)
	{
	   // printf("%lf ",sampleMat.at<float>(j));
	}
	//printf("\n");
	for(int t=0;t<7;t++)
	{
	    sampleMat.at<float>(t+21)=array[34][t+21]*(1.001+x);
	}
	double response=SVM.predict(sampleMat);
	for(int t=0;t<7;t++)
	{
	    sampleMat.at<float>(t+21)=array[34][t+21]*(1+x);
	}
	double response2=SVM.predict(sampleMat);
	double dx=(response2-response)/0.001;
	if(dx>dymax)
	{
	    dymax=dx;
	    ymax=x;
	}
	//printf("%lf\t%lf\t%lf\n",x,dx,response);
    }
    //printf("ymax=%lf,dy=%lf\n",ymax,dymax);
    
    for(int j=0;j<FC;j++)
    {
	sampleMat.at<float>(j)=array[34][j];
    }
    double dzmax,zmax;
    for(double x=0;x>-1;x=x-0.001)
    {
	for(int j=0;j<FC;j++)
	{
	   // printf("%lf ",sampleMat.at<float>(j));
	}
	//printf("\n");
	for(int t=0;t<7;t++)
	{
	    sampleMat.at<float>(t+0)=array[34][t+0]*(1.001+x);
	}
	double response=SVM.predict(sampleMat);
	for(int t=0;t<7;t++)
	{
	    sampleMat.at<float>(t+0)=array[34][t+0]*(1+x);
	}
	double response2=SVM.predict(sampleMat);
	double dx=(response2-response)/0.001;
	if(dx>dzmax)
	{
	    dzmax=dx;
	    zmax=x;
	}
	//printf("%lf\t%lf\t%lf\n",x,dx,response);
    }
    //printf("zmax=%lf,dz=%lf\n",zmax,dzmax); 
    // optimization
    
    real_1d_array x;
    real_1d_array bndl,bndu;
    double x0[FC];
    for (int i=0;i<FC;i++)
    {
        x00[i]=array[34][i];
    }
    for(int i=0;i<7;i++)
	x0[i]=0;
    x0[2]=xmax*3/4;x0[3]=ymax*3/4;
    x0[0]=0.5;x0[1]=0;
    x.setcontent(7,x0);
    double l[7],u[7];
    
    //Smooth		//Sat of H	//Stddev		//frecks	
    l[0]=-0.7;		l[1]=0;		l[2]=x0[2];		l[3]=x0[3];
    u[0]=-0.3;		u[1]=0.5;	u[2]=xmax*2/3; 		u[3]=ymax*2/3;
    
    // Hue		Light	   	Sat
    l[4]=-0.07;		l[5]=-0.12;	l[6]=-0.05;
    u[4]=0.03; 		u[5]=0.05; 	u[6]=0.05;
    
    for(int j=0; j<7; j++)
    {
       rate[j]=x0[j];
	   brate[j]=x0[j];
    }
    
    
    //l[0]=-1;l[1]=0;l[2]=-1;l[3]=-1;
    //u[0]=0;u[0]=1;u[2]=0;u[3]=0;
    
    bndl.setcontent(7,l);
    bndu.setcontent(7,u);
    
    
    double epsg = 0.000000001;
    double epsf = 0;
    double epsx = 0;
    ae_int_t maxits = 0;
    minlmstate state;
    minlmreport rep;

    minlmcreatev(1, x, 0.001, state);
    minlmsetbc(state, bndl, bndu);
    minlmsetcond(state, epsg, epsf, epsx, maxits);
    
    alglib::minlmoptimize(state, function1_fvec);     
    minlmresults(state, x, rep);

    
//     printf("%d\n", int(rep.terminationtype)); // EXPECTED: 4
//     printf("%s\n", x.tostring(6).c_str()); // EXPECTED: [-1,+1]
    
    //printf("[");
//    for(int j=0;j<FC;j++)
// 	printf("%lf ",array[34][j]);
//     printf("]\n");
    
    for(int j=0; j<7; j++)
    {
       double *xx = x.getcontent();
       rate[j]=*(xx+j);
    }
    
    //cerr << rate[0] << " " << rate[1] << " " << rate[2] << " " << rate[3] << " " << rate[4] << " " << rate[5] << " " << rate[6] << endl;
    
    
    for(int t=0;t<7;t++)
	rate[t]=rate[t]*(1 + argv6);
    
    //cerr << rate[0] << " " << rate[1] << " " << rate[2] << " " << rate[3] << " " << rate[4] << " " << rate[5] << " " << rate[6] << endl;
    
    Mat src,src_h,src_m,dest;
	mask=imread(MASK_JPG);
	maskR=imread(MASK3_JPG);
	src_h=imread(imagePath_h);
	src_m=imread(imagePath_m);
	src=imread(imagePath_d);
// 	imshow("test",src_h);
// 	waitKey(0);
// 	imshow("test",src);
// 	waitKey(0);
// 	imshow("test",src_m);
// 	waitKey(0);
	int srate=3;
	resize(src_h,src_h,Size(1536/srate,1536/srate),0,0,CV_INTER_LANCZOS4);
	resize(src,src,Size(1536/srate,1536/srate),0,0,CV_INTER_LANCZOS4);
	resize(src_m,src_m,Size(1536/srate,1536/srate),0,0,CV_INTER_LANCZOS4);
	resize(mask,mask,Size(1536/srate,1536/srate),0,0,CV_INTER_LANCZOS4);
	resize(maskR,maskR,Size(1536/srate,1536/srate),0,0,CV_INTER_LANCZOS4);
	src_h.copyTo(src_hb);
	src_m.copyTo(src_mb);
	src.copyTo(srcb);
	//printf("%d %d\n",src_m.rows,src_m.cols);
    // resynthesis
	clock_t totaltime=clock();
	#pragma omp parallel
	{
		/*while(currEnt>destEnt&&itercount<35)
		{
		GaussianBlur( src_h, dest, Size( 2*itercount+1, 2*itercount+1 ), 0, 0 );
		currEnt=calcEnt(dest,mask);
		itercount+=1;
	   
		}
		printf("currEnt=%lf,K=%d\n",currEnt,2*itercount-1);
		if(itercount!=0)
		GaussianBlur(src_h, src_h, Size( 2*itercount-1, 2*itercount-1 ), 0, 0 );*/
		#pragma omp sections 
		{
			#pragma omp section
			{
				clock_t begin_time = clock();
	
				vector <Mat> mv;
				int size=20*fabs(rate[0]);
				//gpu::GpuMat gpuSrcH(src_h);
				//gpu::GaussianBlur(gpuSrcH,gpuSrcH,Size( 31, 31 ), 0, 0 );
				GaussianBlur(src_h, src_h, Size( 2*size+1, 2*size+1 ), 0, 0 );
				//std::cout << "Elapsed time / GaussianBlur: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
		
				//gpuSrcH.download(src_h);
				begin_time=clock();
				cvtColor(src_h,src_h,CV_BGR2HLS);
				split(src_h,mv);
				mv[2].convertTo(mv[2],-1,rate[1]+1,0);
				merge(mv,src_h);
				cvtColor(src_h,src_h,CV_HLS2BGR);
				//std::cout << "Elapsed time / TestCopy : "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
			}
			#pragma omp section
			{
				clock_t begin_time=clock();
				//printf("===full aFter===\n");
				// full image remove
				{
					Mat mask2(src_m.rows,src_m.cols,CV_8UC1);
					cvtColor(mask,mask2,CV_BGR2GRAY);
					Scalar stdev,mean;
					Mat dest(src_m.rows,src_m.cols,CV_8UC1);
					cvtColor(src_m,dest,CV_BGR2GRAY);
					meanStdDev(dest,mean,stdev,mask2);    
					double pixThres=mean.val[0]-2*stdev.val[0];
					double desThres=mean.val[0]-stdev.val[0];
					double pixMean[3];
					double pixStdev[3];
					meanStdDev(dest,mean,stdev,mask2); 
					//printf("Mean:%lf Stdev:%lf,Thres=%lf\n",mean.val[0],stdev.val[0],pixThres);
	
					float range[] = { 0, (int) pixThres } ;
					int histSize = (int) pixThres;
					const float* histRange = { range };
					bool uniform = true; bool accumulate = false;

					Mat hist;
					meanStdDev(mask2,mean,stdev); 
					double wcount=mean.val[0]/255.0*(src_m.rows*src_m.cols);
	
					calcHist( &dest, 1, 0, mask2, hist, 1, &histSize, &histRange, uniform, accumulate );
					//printf("%d\n",hist.rows);
					double aveG=0;
					int xt=0;
					double tG=0;
					for(int q=0;q<hist.rows;q++)
						tG+=hist.at<float>(q)/wcount;
	    
					for(int q=0;q<hist.rows;q++)
					{
						aveG+=hist.at<float>(q)/wcount;
						//printf("q:%d,%lf\n",q,aveG);
						if(aveG>=(1-fabs(rate[3]))*tG)
						{
							//printf("DESTRATE:%lf,q=%d\n",1-fabs(rate[3]),q);
							xt=q;
							break;
						}
		    
					}
					//std::cout << "Elapsed time / Histogram: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
    
					begin_time=clock();
					//Mat dest2;
					//dest.copyTo(dest2);
					#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						uchar *t;
						Vec3b *s;
						t=dest.ptr<uchar>(x);
						s=mask2.ptr<Vec3b>(x);
						Vec3b *w;
						w=src_m.ptr<Vec3b>(x);
						for(int y=0;y<src_m.cols;y++)
						{
							//printf("%d %d\n",x,y);
							if(t[y]<xt && s[y][0]!=0)
							{
								//t[y]=t[y]*pixThres/xt;
								for(int v=0;v<3;v++)
									w[y][v]=w[y][v]*pixThres/xt;
								//printf("Yeah: %lf\n",(double)pixThres/(double)xt);
							}
							if(t[y]>=xt && t[y]<pixThres && s[y][0]!=0)
							{
								//t[y]=pixThres;
								double dg=0.2126*w[y][0]+0.7152*w[y][1]+0.0722*w[y][2];
								for(int v=0;v<3;v++)
									w[y][v]=w[y][v]*pixThres/dg;
							}
						}
					}
					//getchar();
					/*#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						for(int y=0;y<src_m.cols;y++)
						{
							for(int t=0;t<3;t++)
							{
								if(dest2.at<uchar>(x,y)!=0)
									src_m.at<Vec3b>(x,y)[t]=src_m.at<Vec3b>(x,y)[t]*dest.at<uchar>(x,y)/dest2.at<uchar>(x,y);
								else
									src_m.at<Vec3b>(x,y)[t]=0;
							}
						}
					}*/
					//std::cout << "Elapsed time / frecks: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
    
					begin_time=clock();
					//dest.copyTo(dest2);
					meanStdDev(dest,mean,stdev,mask2); 
					double meanval=mean.val[0];
					#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						Vec3b *w;
						w=src_m.ptr<Vec3b>(x);
						uchar *v;
						v=dest.ptr<uchar>(x);
						for(int y=0;y<src_m.cols;y++)
						{
							double t=meanval-(meanval-v[y])*(1+rate[2]);
							//double dg=0.2126*w[y][0]+0.7152*w[y][1]+0.0722*w[y][2];
							for(int j=0;j<3;j++)
								w[y][j]=w[y][j]*t/v[y];

						}
					}
					/*#pragma omp parallel for num_threads(8)
					for(int x=0;x<src_m.rows;x++)
					{
						for(int y=0;y<src_m.cols;y++)
						{
							for(int t=0;t<3;t++)
							{
								if(dest2.at<uchar>(x,y)!=0)
									src_m.at<Vec3b>(x,y)[t]=src_m.at<Vec3b>(x,y)[t]*dest.at<uchar>(x,y)/dest2.at<uchar>(x,y);
								else
									src_m.at<Vec3b>(x,y)[t]=0;
							}
						}
					}*/
					//std::cout << "Elapsed time / stddev: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
				}
			}
		}
	 }
		//imshow("src_m",src_m);
		//waitKey(0);
		clock_t begin_time=clock();
		//printf("Maybe Here?\n");
		Mat dif=src_h;
		Mat ddif=src_m;
		Mat dst(dif.rows,dif.cols,CV_8UC3);
		double smax=393.088;
		double smin=-0.969;
		double dmax=97.104;
		double dmin=-76.391;
		double ddmax=90;
		double ddmin=-90;
		/*Mat doubleh(src_h.rows,src_h.cols,CV_32FC3);
		Mat doublem(src_h.rows,src_h.cols,CV_32FC3);
		Mat doubled(src_h.rows,src_h.cols,CV_32FC3);
		clock_t testtime=clock();
		#pragma omp parallel 
		{
			#pragma omp sections
			{
				#pragma omp section
				{
					clock_t sectime=clock();
					src_h.convertTo(doubleh,CV_32FC3);
					doubleh.convertTo(doubleh,-1,(dmax-dmin)/255,dmin);
					std::cout << "Elapsed time / total: "  << float( clock () - sectime ) /  CLOCKS_PER_SEC << endl;	
				}
				#pragma omp section
				{
					clock_t sectime=clock();
					src_m.convertTo(doublem,CV_32FC3);
					doublem.convertTo(doublem,-1,(ddmax-ddmin)/255,ddmin);
					std::cout << "Elapsed time / total: "  << float( clock () - sectime ) /  CLOCKS_PER_SEC << endl;	
				}
				#pragma omp section
				{
					clock_t sectime=clock();
					src.convertTo(doubled,CV_32FC3);
					doubled.convertTo(doubled,-1,(smax-smin)/255,smin);
					std::cout << "Elapsed time / total: "  << float( clock () - sectime ) /  CLOCKS_PER_SEC << endl;	
				}
			}
		}
		std::cout << "Elapsed time / total: "  << float( clock () - testtime ) /  CLOCKS_PER_SEC << endl;
		clock_t sectime=clock();
		add(doubleh,doublem,doubleh,Mat(),-1);
		std::cout << "Elapsed time / ADD: "  << float( clock () - testtime ) /  CLOCKS_PER_SEC << endl;
		add(doubleh,doubled,doubled,Mat(),-1);
			doubled.copyTo(dst);
		*/
		Vec3b *p,*q,*a,*b;
		#pragma omp parallel for num_threads(8) private(p,q,a,b)
		for(int i=0;i<src.rows;i++)
		{
		  for(int j=0;j<src.cols;j++)
		  {
			  p = src.ptr<Vec3b>(i);
			  q = dif.ptr<Vec3b>(i);
			  a = ddif.ptr<Vec3b>(i);
			  b = dst.ptr<Vec3b>(i);
			  for(int t=0;t<3;t++)
			  {
				double l=p[j][t];
				double r=q[j][t];
				double m=a[j][t];
				l=l/255*(smax-smin)+smin;
				r=r/255*(dmax-dmin)+dmin;
				m=m/255*(ddmax-ddmin)+ddmin;
				if(l+r+m>255)
					b[j][t]=255;
				else
				{
					if(l+r+m<0)
						b[j][t]=0;
					else
						b[j][t]=l+r+m;
				}
			  }
		   }
		}
		//std::cout << "Elapsed time / trimerge: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;
    
    //system("./trimerge ./r34_d.jpg ./final_h.jpg ./final_m.jpg 393.088 -0.969 97.104 -76.391 90 -90");
	begin_time=clock();
    Mat corrSrc=dst;
    cvtColor(corrSrc,corrSrc,CV_BGR2HLS);
	/*vector <Mat> mv;
	split(corrSrc,mv);
	#pragma omp parallel for
	for(int i=0;i<3;i++)
		mv[i].convertTo(mv[i],-1,1+rate[4+i],0);
	merge(mv,corrSrc);*/
	
	#pragma omp parallel for num_threads(8) private(p,q,a,b)
    for(int x=0;x<corrSrc.rows;x++)
    {
		p = corrSrc.ptr<Vec3b>(x);
		q = src.ptr<Vec3b>(x);
		for(int y=0;y<corrSrc.cols;y++)
		{
			if(maskR.at<Vec3b>(x,y)[0])
			{
				for(int j=0;j<3;j++)
				{
					int r=p[y][j];
					if(r*(1+rate[j+4])<255)
						p[y][j]=r*(1+rate[j+4]);
					else
						p[y][j]=255;
				}
			}
		}
    }
    cvtColor(corrSrc,corrSrc,CV_HLS2BGR);

	return corrSrc;

	//std::cout << "Elapsed time / other: "  << float( clock () - begin_time ) /  CLOCKS_PER_SEC << endl;	
	//std::cout << "Elapsed time / total: "  << float( clock () - totaltime ) /  CLOCKS_PER_SEC << endl;	
    //corrSrc.convertTo(corrSrc,-1,0.93,15);
// 	namedWindow("Result",CV_WINDOW_AUTOSIZE);
// 	createTrackbar("rate", "Result", 0, 100, scrollbar_onchange, &rate);
// 	imshow("Result",corrSrc);
// 	waitKey(0);
	
// 	for(int j=0;j<100;j++)
// 	{
// 		//printf("====j=%d====\n",j);
// 		scrollbar_onchange(j,0);
// 		waitKey(50);
// 	}
// 
//     imwrite("FINAL_F2.jpg",corrSrc);
//     
//     //printf("==============prediction===============");
//     for(int j=0;j<FC;j++)
//     {
// 	sampleMat.at<float>(j)=array[34][j];
//     }
//     for (int j=0;j<7;j++)
//     {
// 	sampleMat.at<float>(j)=sampleMat.at<float>(j)*(1+rate[0]);
//     }
//      for (int j=7;j<14;j++)
//     {
// 	sampleMat.at<float>(j)=sampleMat.at<float>(j)*(1+rate[1]);
//     }
//      for (int j=14;j<21;j++)
//     {
// 	sampleMat.at<float>(j)=sampleMat.at<float>(j)*(1+rate[2]);
//     }
//      for (int j=21;j<28;j++)
//     {
// 	sampleMat.at<float>(j)=sampleMat.at<float>(j)*(1+rate[3]);
//     }
//     sampleMat.at<float>(28)=sampleMat.at<float>(28)*(1+rate[4]);
//     sampleMat.at<float>(29)=sampleMat.at<float>(29)*(1+rate[5]);
//     sampleMat.at<float>(30)=sampleMat.at<float>(30)*(1+rate[6]);
//     
//     //printf("\n%lf\n",SVM.predict(sampleMat));

//    return 0;
}

int main(){
	Mat mat;
	mat = main_process(0.5, 1, IMAGE_H, IMAGE_M, IMAGE_D);
	imshow("result", mat);
	cvWaitKey(0);
	return 1;
}


